
function [p q] = implicit_theta(dt,beta)
global U  h  g  p  q  nt dx x h1 p_global;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;

%---- Construct the spatial Discretization Matrix------%
A = zeros(x-1,x-1);
rhs = zeros(2*x-1,1);

%A(1,1) = 1;
%A(1,x) = -U*0.5*dx1; %%% Not periodic Boundary condition but Dirichlet
A(1,2) =  U*0.5*dx1;
for k=2:x-2
    A(k,k-1) = -U*0.5*dx1;
    A(k,k+1) =  U*0.5*dx1;
end

A(x-1,x-2) = -U*dx1;
A(x-1,x-1)   =  U*dx1;

L=zeros(2*x-1,2*x-1);
L(1:x-1,1:x-1)=A;

B = zeros(x,x);

%Neumann Boundary Conditions
B(1,1) = -U*0.5*dx1;
B(1,2) =  U*0.5*dx1;
for k=2:x-1
    B(k,k-1) = -U*0.5*dx1;
    B(k,k+1) =  U*0.5*dx1;
end

B(x,x-1) = -U*0.5*dx1;
B(x,x)   =  U*0.5*dx1;

L(x:2*x-1,x:2*x-1)=B;

A = zeros(x-1,x);

for k=1:x-2
    A(k,k)= h*dx2;
    A(k,k+1)= -2*h*dx2;
    A(k,k+2) = h*dx2;
end

A(x-1,x-1)= h*dx2;
A(x-1,x) = -h*dx2;

L(1:x-1,x:2*x-1)=A;

for k=1:x-1
   L(x+k,k)=g;
end 

I = eye(2*x-1);
% Periodic Boundary Conditions
sol=rhs;
sol(1:x-1,1) =p;
sol(x:2*x-1,1) =q;
%p_global(:,1)=p;

for n=2:nt+1;
    gamma_bc =1.0;%sin(2*nt*dt);
    rhs = (dt2*I - (1-beta)*L)*sol;
    rhs(1) = rhs(1) + U*0.5*dx1*gamma_bc;
    rhs(x) = rhs(x) - g*gamma_bc;
    sol=(dt2*I+beta*L)\rhs;    
    p = sol(1:x-1,1);
    q= sol(x:2*x-1,1);
%p_global(:,n)=p;
%     if rem(n,100)==0
%         
%         %refreshdata(h1,'caller') % Evaluate p in the function workspace
%      %   drawnow
%     end
    
end
display('Completed Successfully');